Early dust evolution in protostellar accretion disks 
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^ : ABSTRACT 

\ We investigate dust dynamics and evolution during the formation of a protostellar accretion 

'sj" ' disk around intermediate mass stars via 2D numerical simulations. Using three different detailed 

dust models, compact spherical particles, fractal BPCA grains, and BCCA grains, we find that 
even during the early collapse and the first ~ 10** yr of dynamical disk evolution, the initial dust 
\ size distribution is strongly modified. Close to the disk's midplane coagulation produces dust 

particles of sizes of several 10 /im (for compact spherical grains) up to several mm (for fluffy 
BCCA grains), whereas in the vicinity of the accretion shock front (located several density scale 
heights above the disk), large velocity differences inhibit coagulation. Dust particles larger than 
O \ about 1 /Ltm segregate from the smaller grains behind the accretion shock. Due to the combined 

4-J ' effects of coagulation and grain segregation the infrared dust emission is modified. Throughout 

\ the accretion disk a MRN dust distribution provides a poor description of the general dust 

properties. Estimates of the consequences of the "freezing out" of molecules in protostellar 
disks should consider strongly modified grains. Physical model parameters such as the limiting 
' sticking strength and the grains' resistivity against shattering are crucial factors determining the 

^ \ degree of coagulation reached. In dense regions (e.g. in the mid-plane of the disk) a steady- 

state is quickly attained; for the parameters used here the coagulation time scale for 0.1 /im dust 
particles is ~ 1 yr (10~^^g cm~^ / g). High above the equatorial plane coagulation equilibrium is 
not reached due to the much lower densities. Here, the dust size distribution is affected primarily 
by differential advection, rather than coagulation. The influence of grain evolution and grain 
dynamics on the disk's near infrared continuum appearance during the disk's formation phase is 
only slight, because the most strongly coagulated grains are embedded deep within the accretion 
disk. 

Subject headings: accretion, accretion disks — hydrodynamics — radiative transfer — solar system: 
formation — stars: formation 

1. Introduction as their observable appearance. Dust provides the 

seeds for planetesimals, which in turn evolve into 
Dust in protostellar envelopes and accretion constituents of a planetary system: comets, 

disks is a major component of pre-stellar matter, planets, moons, and the debris associated with as- 

strongly influencing the thermodynamical and gas Ceroids and Kuiper belt objects. Interstellar dust 

dynamical behavior of these young objects as well evolves significantly from its initial state in re- 
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cently formed molecular clouds up to the forma- 
tion of planetesimals around stars and, as the gas 
density increases, it does so at an ever increasing 
rate. In molecular clouds it takes several 10^ yr 
to build up large fluffy grains by dust coagulation 
(Ossenkopf 1993; Weidenschilling & Ruzmaikina 
1994) . In the midplane of accretion disks this time 
scale shrinks to about 10^ yr due to the high den- 
sities there (Mizuno, Markiewicz, & Volk 1988; 
Mizuno 1989; Schmitt, Henning, & Mucha 1997). 

It would be naive to imply, however, that it is 
only a matter of time before coagulation produces 
the first planetesimals of several km in size. Other 
processes affect grain growth and evolution: or- 
bital decay, shattering, cratering, sputtering, and 
compacting of amorphous grains, in addition to 
adsorption, outgassing, and chemical reprocessing 
of molecules. These processes depend critically 
on the physical conditions within the disk and on 
the grains' relative velocities. Large compact dust 
grains ( > 10 /im) can decouple from the gas and 
gain large relative velocities to each other which 
could prevent coagulation due to a limited stick- 
ing strength. By contrast, large fractal dust grains 
— such as Ballistic Cluster Cluster Agglomeration 
(BCCA) particles — are always well coupled to 
the gas component; large BCCA particles do not 
achieve sufficient relative velocities to coagulate 
effectively (Schmitt et al. 1997). The assumption 
of an "intermediate" type of fractal grains, the 
so called Ballistic Particle Cluster Agglomeration 
(BPCA) particles, could alleviate this problem. 
Such coagulates possibly couple sufficiently well to 
the gas to prevent high relative velocities and be- 
have in the limit of a large number of constituent 
grains like compact particles, so that turbulence 
and systematic relative velocities can drive coag- 
ulation (Ossenkopf 1993). 

A high sticking strength is a necessary pre- 
requisite for effective coagulation up to planetes- 
imal sizes, since turbulent relative velocities can 
prevent the coagulation of pre-planetesimal dust 
grains (Weidenschilling & Cuzzi 1993) . This could 
be achieved with ice or frost layers on the colliding 
particles' surfaces (Bridges et al. 1996; Supulver 
et al. 1997). 

Dust coagulation leads to important modifica- 
tions of the protostellar matter. Turbulence in 
accretion disks is strongly coupled to the opacity 
of the medium which in turn is dependent on the 



type of dust material and the dust grain size distri- 
bution. A high degree of coagulation implies sig- 
nificant reduction of the dust opacity (Mizuno et 
al. 1988). Reduced opacity is necessary to damp 
turbulence which otherwise would be a mayor ob- 
stacle for the pre-planetesimal dust particles to 
settle down to the disk's equatorial plane (Weiden- 
schilling 1984). During the formation of an accre- 
tion disk, however, virtually unprocessed dust ma- 
terial is continuously being supplied by the parent 
molecular cloud core. "Second generation" coagu- 
lation occurs (Mizuno 1989) which also influences 
opacity and the dynamics of the disk. 

For this study we calculated the dust evolution 
during the formation of a protostellar accretion 
disk using a multi-component radiation hydrody- 
namics code, an improved version of a RHD-code 
designed for one component (Yorke, Bodenheimer, 
& Laughlin 1995; Yorke & Bodenheimer 1999). 
Different dust models are applied and the influ- 
ence of the dust-dust sticking strength is investi- 
gated. Brownian motion, turbulence, differential 
radiative acceleration and gravitative sedimenta- 
tion with size dependent relaxation time scales are 
considered as sources of relative velocities between 
the dust grains. Mean dust opacities are calcu- 
lated directly from the actual dust size distribu- 
tion and are continuously updated in the radiation 
transport module. Finally, a diagnostic frequency 
and angle dependent radiation transport code is 
used to produce synthetic dust continuum emis- 
sion maps and spectra which give information on 
observational consequences. 

In section 2 we describe our models for the dust 
and in section 3 we sketch our mimca'ical approach 
to the problem. The initial conditions for the par- 
ticiilar cases considered are introduced in section 
4 and the numerical results of the simulations are 
presented in section 5. In the final section 6 these 
results are discussed in light of observable conse- 
quences and conclusions are drawn. 

2. Physics of dust grains in star forming 
regions 

According to our present picture of the inter- 
stellar medium, dust grains more or less uniformly 
permeate gas clouds, contributing about one per- 
cent of their total mass. The grain size distri- 
bution follows a power law with an exponent of 
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about -3.5 (Mathis, Rumpl, & Nordskx'k 1977; 
hereafter MRN). MRN estimated maximum and 
minimum grain sizes at several nm and several /xm, 
respectively. The constituent dust grains wore as- 
sumed to be mainly "astrophysical" silicates with 
some contribution by "graphites" . With this dust 
model they could fit the overall extinction curve 
of the diffuse interstellar medium quite well. How- 
ever, observations of dense molecular cloud cores 
show that the upper limiting grain size is shifted to 
larger radii (polarization measurements of Vrba, 
Coyne, & Tapia 1993; theoretical predictions by 
Fischer, Hcnning, & Yorke 1994). Coagulation is 
assumed to grow large grains in this environment. 

Three different dust models are considered. 
First, we treat the dust as compact spherical par- 
ticles. For this model there are a several the- 
oretical studies which deal with dust opacities 
(Yorke 1988), sticking strengths (Chokshi, Tie- 
lens, & Hollenbach 1993), the physics of dust shat- 
tering (Jones, Tielens, & Hollenbach 1996), and 
gas-dust coupling strengths (Yorke 1979). Thus, 
the basic physical properties are reasonably well 
defined, although more recent experimental work 
indicates that the critical sticking velocities tend 
to be higher than theoretical estimates by about 
a factor of 10 (Poppe & Blum 1997). 

Wc next consider fractal BPCA grains and 
BCCA grains. Osscnkopf (1993) developed a the- 
oretical dust model for these fluffy particles and 
provided analytical expressions for the gas-dust 
and the dust-dust interaction cross sections as a 
function of the compact volume of the grains. 
Henning & Stognienko (1996) have calculated 
opacities for those fractal dust particles and Wurm 
(1997) generalized the sticking strength of Chok- 
shi et al. (1993) to fractal coagulates. Again, 
recent experimental studies indicate critical stick- 
ing velocities to be about an order of magnitude 
larger than the theoretical values (Poppe & Blum 
1997). 

2.1. Compact spherical grains 

2.1.1. Opacities and radiative acceleration 

Because dust is the predominant source of ex- 
tinction in the temperature, density, and wave- 
length regimes considered (e.g. Yorke &: Henning 
1994), we have included only the dust's contribu- 
tion to opacity and radiative acceleration. The 




Fig. 1. — Specific extinction coefScient of com- 
pact spherical dust grains with sizes of 5 nm (solid 
line), 0.1 /xm {dotted line), bfim {dashed line) and 
0.2mm {dot-dashed line). 

specific extinction K^j of the component i at wave- 
length A is proportional to the specific cross sec- 
tion of the spherical grains: 

'^ti = Ql>o,Vm, (1) 

JV JV 

= (2) 

i=l i=l 
<,i = «0 (l-^A,i5A,i) (3) 

The extinction gain factor Q^xl^ the albedo A\^i, 
and the asymmetry parameter gx^i = (cos 8 a) j 
are calculated from Mie-theory using the dielec- 
tric constants for "astronomical" silicates given 
by Draine & Lee (1984). The net specific extinc- 
tion of the dust grains k'^*' is the weighted sum 
over the different constituent grain sizes (Fig. 1). 
Here, ^ is the specific radiation pressure opac- 
ity needed to determine the interaction of the dust 
particles with the stellar radiation flux Y\. The 
radiative acceleration of the dust particles is given 

2.1.2. Sticking of grains 

Sticking of the spherical dust grains can be de- 
scribed by application of the theory of elasticity 
(Chokshi et al. 1993). One determines the bind- 
ing energy of two spheres brought together under 
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consideration of clastic wave dissipation in the two 
bodies. The resultant binding energy is set equal 
to the relative kinetic energy of the two colliding 
particles. This gives a critical sticking velocity, 
above which coagulation ceases and the particles 
bounce. Thus, t^stick is given by: 

= 9.6 7'/' R*/' E-V3 (5) 
Vstick = \/2£;stick/M (6) 

The constants 7 and E denote the surface en- 
ergy per unit and Young's modulus of the 
dust material, respectively. R = rir2/(ri + r2) 
is the reduced radius of the contact surfaces and 
H = mim2/ (toi + 1712) is the reduced mass of the 
two colliding dust particles. This expression differs 
by a numerical factor of about 3 from the original 
formula by Chokshi et al. (1993) when applied to 
spheres, because a correction of Dominik & Tielens 
(1997) has been applied (see discussion by Wurm 
1997). Note that Beckwith, Henning, and Nak- 
agawa (2000) give a slightly different formula for 
^^stick (quoting the same authors as above): 



^^stick = 1-07 



7' 



5/6 



R5/6 El/3 



(7) 



where pg is the specific mass of the grain material. 
The formulae 5/6 and Eq. 7 arc identical when 
ri = r2 and mi = m2 but differ by 18% (50%) 
when the mass ratio of the colliding particles is 10 
(100) and by a factor of almost two for extreme 
dust mass ratios. 

An ice layer on the grains' surfaces enhances 
the critical sticking velocity by more than an order 
of magnitude compared to a pure silicate surface. 
The ice surfaces arc destroyed at dust tempera- 
tures of about 125 K. Thus, pure silicate grains 
are assumed where; the dust temperature exceeds 
the ice melting limit. According to the experi- 
ments of Poppe & Blum (1997) the critical stick- 
ing velocities are increased by a factor of 10 for 
the simulations. 

2.1.3. Grain shattering 

When the relative velocity rises above a criti- 
cal velocity Vcnt, the dust particles are shattered. 
This is a gradual process, starting with crater for- 
mation on the grains' surfaces and ending with 
total disruption of projectile and target. We as- 
sume a critical shattering velocity Vcrit for silicates 



(see Jones et al. 1996): 

t^crit = 2.7 km s~^ 

Jones et al. (1996) also give analytic expres- 
sions for the ejected mass during a shattering en- 
counter. When half of the mass of the target par- 
ticle (the larger of two colliding dust grains) is 
shocked, they assume total disruption; the debris 
particles are assumed to follow a power law size 
distribution with an exponent of —3.0 to —3.5. 
In our simulations we choose an exponent —3.5. 
The upper size limit of the fragments first grows 
with increasing collision velocity (i.e. increasing 
crater volume) and starts to decrease inversely 
proportional to the relative velocity when disrup- 
tion starts (see appendix). 

2.1.4. Gas dust interaction 

The dust grains are coupled to the gas by dy- 
namical friction. According to the gas densities in 
star forming regions the Epstein coupling law (Ep- 
stein 1923; Weidenschilling 1977) is valid because 
the mean free path A/ of the gas molecules is large 
compared to the radii of the dust particles: 



A/ 



^ /^gas^gas 

10'" cm 



10-"gcm-3 



>a, (8) 



Here, we use an extension of the Epstein law for su- 
perthermal relative velocities (Yorke 1979), which 
often occur in later stages of the protostellar col- 
lapse: 



'dt 



^e^V^TW^iv-v,) (9) 



The interaction term (9vi/9f)ww describes the 
coupling of dust particles with cross section (Tj 
and mass mj with gas of density g and isother- 
mal sound speed Cg. 

2.2. Fractal grains 

2.2.1. Mass to radius relation 

There is no simple relation between mass and 
radius for fractal grains as there is for compact 
spherical dust particles. Here, a grain model which 
describes the transition from the compact con- 
stituent grains (at the lower radius limit of the size 
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distribution) to the fractals (at the higher end) 
must be appHed. For BPCA and BCCA grains 
Henning & Stognienko (1993) give an analytic ex- 
pression for the filling factor / in relation to the 
extremal radius rex,i of the grains: 



An 



3 

— ^'bulk fex,i 



fi 



fP<^^ = 0.0457 



jCCA ^ Q 279 



1 + 696 



0.01 /xm 
1 + 4.01 



0.01 ;um 

-1.04 



-3.93 



0.01 Aim 



-1.34' 



(10) 



(11) 



(12) 



The extremal radius rex,i is defined as the radius 
of the minimal envelope sphere covering a fractal 
particle with bulk density ^pbuik- The constituent 
grains are assumed to have radii of 0.01 fim. 

2.2.2. Opacities and radiative acceleration 




Fig. 2. — Specific extinction coefficient of fractal 
dust particles with sizes of 5 nm [solid line), 0.1 
{dotted line), 5/Ltm {dashed line), and 0.2mm {dot- 
dashed line). 

Dust opacities for coagulated grains, averaged 
over a given size distribution, have been calculated 
by Henning & Stognienko (1996). R. Schrapler 
(private communication), with the authors' per- 
mission kindly provided us with their basic size 
dependent specific extinction coefficients. Wc used 
the opacities for olivine ([Fe,Mg]2Si04) dust grains 
(Fig. 2). 



2.2.3. Sticking of grains 

In contrast to compact spherical grains the frac- 
tal coagulates stick at their limb structures which 
are built by the small constituent grains. The re- 
duced radius 'R' in Eq. 5 for the critical sticking 
energy must be calculated using the radii of these 
constituent grains. In this manner the formula 
for the sticking velocity by Chokshi ct al. (1993) 
can be generalized to fractal particles (see section 
2.1.2). 

2.2.4- Grain shattering 

Two approaches were followed. First, grain 
shattering is treated the same way as for the com- 
pact spherical grains. The analysis depends only 
on material parameters, the relative velocity, and 
the masses of the colliding particles. Although 
this assumption may seem inadequate, it poses 
an upper limit to the grains' resistance against 
shattering. Fractals seem to be rather fragile con- 
structs with low binding energy when only van 
der Waals adhesion is considered. However, ex- 
perimental studies show that fractal particles are 
always more resistant against destruction than the 
predictions of theoretical models (Poppe & Blum 
1997; Wurm 1997). Fractals possess a multitude 
of vibrational excitation modes which provide a 
wide channel to dissipate the kinetic energy of the 
impacting grain. 

As a second approach we adapt the shattering 
model of Dominik & Tielens (1997): The criti- 
cal shattering velocity Vait is proportional to the 
sticking velocity VgUck and the number of contact 
points between the two colliding grains. The val- 
ues thus obtained are rather low compared to the 
model of Jones et al. (1996). 

2.2.5. Gas-dust interaction 

The gas-dust interaction also has to be modi- 
fied due to the fractal structure of BPCA grains. 
In order to find an analytic expression for the effec- 
tive cross section, Ossenkopf (1993) fitted numeri- 
cal calculations of the cross section a of coagulates 
in dependence of their compact volume V: 



a 





b 







(13) 



15.2, 
.4.7, 



2 < V/Vq < 1000 
V/Vq > 1000 
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2.86, 2 < y/Fo < 1000 

9.02, y/^o > 1000 

0.096, 2 < y/Vo < 1000 

0.503, y/Fo > 1000 



For BCCA particles the following expression 
has been derived: 



like BPCA 



0.692 (^)° "(l + ^) 



^<25 



(14) 



Here, the normalization factors cto and Vq are 
the cross section and the volume of the (com- 
pact) constituent grains ao = n ■ (0.01 fim)^ and 
Vc) = 47r/3 • (0.01 yum)^. The compact volume 
can be calculated using the filling factor of sec- 
tion 2.2.1. Thus, a relation between Tex and the 
cross section a of the grains can be established and 
used in the interaction term of section 2.1.4. 

2.2.6. Dust-dust interaction 

For the dust-dust interaction (coagulation) yet 
another radius definition has been introduced. 
The toothing radius rtooth is defined as half the av- 
erage distance of the center of mass of two sticking 
grains. It measures the penetration of two frac- 
tal grains and has been determined by Ossenkopf 
(1993) by fits to numerical calculations which sim- 
ulated the scan of the "coastline" of the larger dust 
particle by the smaller one. Thus, the collisional 
cross section for grain-grain collisions acoU reads: 

CTcoll = 1^'^ (1 - C') 



+ [4n-8.3{l-C'-'')]{rilLr (15) 



with: 



tooth/ ' tooth 

r,„„,, = 0.72 y V3 ^0.21 , 1 



0.216 



X 



(16) 
(17) 
(18) 



The superscripts (1) and (2) denote the grain with 
the larger and the smaller toothing radius, respec- 
tively. 

2.3. Coagulation and shattering 

Adding the different source terms of grain ac- 
celeration, it becomes obvious that the dust par- 
ticles will gain relative velocities to each other. 



The absolute values of these velocities depend on 
the dust absorption and scattering cross sections, 
the radiative flux, the gas density and the specific 
cross section of the grains: 
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Vi)2(v-Vi) (19) 



In addition to these systematic relative velocities 
a random contribution is caused by Brownian mo- 
tion and turbulence. Because our hydrodynam- 
ical grid is too coarse to resolve typical turbu- 
lent length scales, and because turbulence needs 
a three dimensional treatment, we apply a tur- 
bulence model developed by Volk et al. (1980). 
Is is coupled to the turbulent angular momentum 
transport (Shakura & Sunyaev 1973) by the pa- 
rameter a. Here, the macroscopic turbulent ve- 
locity I'turb to a fraction a of the isothermal 
sound speed Cg, and the macroscopic revolution 
time scale f^^.^, is proportional to the orbital pe- 
riod f2: 



nurb 



t 



turb 



(20) 
(21) 



A Kolmogorov-type turbulence spectrum is as- 
sumed, whereby the turbulent energy is trans- 
ported from the largest eddies down to the small- 
est at a constant rate. The corresponding scales of 
the smallest eddies are defined where the turbulent 
Reynolds number Rcg of the gas with viscosity r] 
equals one (i.e. the flow becomes laminar, Lang 
1974): 



''turb 



""turb 



"turb -^■'^0 



""turb 
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S ^turb '^turb 

■n 



turb/ 



''turb 
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(22) 

(23) 

(24) 



The back reaction of the gas turbulence on the 
dust particles depends on the coupling strength 
between grains and gas. This strength can be mea- 
sured by the correlation time scale Tj: 



Tj = |v - Wi\/{dvi/dt)^ 



(25) 
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According to the analytical fit of Wcidcn- 
schilling (1984) to the numerical results of Volk et 
al. (1980), the random relative velocities between 
two grains with correlation time scales ri and T2 
(ti < T2) can be expressed as follows: 



6v, 



turb 



'^turb 



if Tl,T2 < tl^,^ 
ifTi<Crb<^2 



,,0 . tLrhi-ri+T2) 



^turb 

,0 



V, 



2tiT2 

3t2 



turb T1+T2 



ifCb<n. 
otherwise 



(26) 



T2 



The cutoff of the turbulent eddies at the lower size 
end of the eddy spectrum is included according to 
the considerations of Weidenschilling (1984). The 
contribution of Brownian motion to the random 
part of the relative velocities between dust grains 
is only important for small (a, < 1 /xm) grains: 



Sv] 



l8kT', 



■ mo- 



brown 



rriimj 



(27) 



Thus, th{^ dust particles achieve a total rela- 
tive velocity diH,j = (Sv^y^^ + 6v^^^^^ + Svl^^^J^^^, 
where Svgyst denotes the systematic relative veloc- 
ities. These relative velocities are then evaluated 
according to the considerations of section 2.1 or 
2.2. 

If the velocities are sufficiently low, the particles 
can coagulate. This is mathematically described 
by the coagulation equation: 

'dn{m)^ 



dt 



2J J c>:(to,'7i )n{m)n{m ) 
5{m — m — m ) dm dm 



n{m) I a{m, m ) n{m ) dm (28) 



with: 



a{m ,m ) = p acoui'm , w ) Sv{m , m ) (29) 

The variables n(m), icoiii'ni , m ) and 5v{m , m, ) 
are the number density, the relative interaction 
cross section, and the relative velocity of grains 
with masses m, m and m , respectively. The 
sticking probability p controls the onset of bounc- 
ing when the relative velocities become too high 
(p ^ 0). Grain shattering is described by a gen- 
eralization of the above coagulation equation: 



' dn{my 

dt 



shat 



m ) n(m ) n{m ) 



7(m, m ,m , 6v{m , m )) dm dm 
— n{m) J P{m, m ) n(m ) dm (30) 



with: 



l3{m ,m ) = q acoii{m , m ) Sv{m , m ) (31) 

Here, the function 7 distributes the shattered dust 
fragments to the appropriate mass bins (see ap- 
pendix). Again, the shattering probability q con- 
trols the onset of shattering above the critical ve- 
locity (g — > 1). The total dust evolution operator 
(coagulation/shattering) is the sum of both partial 
operators. 

3. Numerical techniques 

To simulate the collapse of a gravitationally 
unstable rotating molecular cloud core we apply 
a multicomponent radiation hydrodynamics code 
with detailed dust dynamics (grain drift, coagula- 
tion, shattering). To keep the problem tractable 
axial symmetry is assumed ("2.5 D" in cylindrical 
coordinates). An explicit nested grid technique is 
applied to resolve the inner parts of the accretion 
disk (Yorke & Kaisig 1995). 

3.1. Solution of the coagulation/shattering 
equation 

To solve the combined coagulation/shattering 
equation numerically the dust size distribution 
is binned into N discrete logarithmically spaced 
mass intervals. The continuous equation (section 
2.3) is therefore discretized: 



dt 



dt 

N N 



drik 
dt 



shat 



i=i 3=1 

N 



(32) 



with: 



1 otherwise 

= ("ifc + TOfe_i)/2 = rrij^^i 
TOj -|- mj 

Qijk ~ 

ruk 

Gk{mi,m,j,Svij) € [0, 1] 



■Gk{mi,mj,6vij) 



(33) 
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N 



^ Gk {rrii , nij , 6vij ) = 1 



k=l 



The distribution of the shattered fragments 
is given by the discrete distribution function 
Gk{mi,mj,6vij) (see appendix). The kernels 
aij and f3ij are the discretized counterparts of 
a{m ,'m ) and /3(m , m ) (section 2.3). Substitut- 
ing backward time differences for ah time deriva- 
tives, this nonlinear integro-differential-equation 
can be brought into the form: 



(n* - n)/Ai = ^(n*)n* 



(34) 



where n = (ni, n2, . . . , n^v) are the (known) den- 
sities at the beginning of the time step and n* are 
the corresponding values after time At which are 
to be determined, ^(n* ) is a matrix which is con- 
structed from the right hand side of Eq. 32. The 
advantages of backward time derivatives are: a) 
the difference equations are numerically stable for 
all choices of time steps At, and b) the numerical 
solution approaches the correct steady-state solu- 
tion as At — > c». We iteratively solve the implicit 
equation 34 for each cell in the numerical grid 
using a multidimensional Newton- Raphson algo- 
rithm. In order to optimize convergence we adap- 
tively reduce the time step with respect to the rel- 
atively large time step used for the explicit hydro- 
dynamics. We tested our solver by calculating the 
solution of simple coagulation problems for which 
analytic solutions exist (Wetherill 1990). The cor- 
respondence was very good at high resolutions of 
mass binning. At resolutions comparable to those 
used in the collapse calculations the accuracy suf- 
fers. Sharply peaked mass distributions and dis- 
continuous mass distributions become more diffuse 
with the passage of time. We do not consider this 
to be a serious flaw, however, because of the large 
uncertainties of the details of grain growth and 
destruction. The total mass of the dust compo- 
nent was always conserved within rounding errors 
(Suttner, Yorke, & Lin 1999). 

3.2. The multicomponent radiation hydro- 
dyneimic code 

An explicit /implicit method is used to solve the 
coupled hydrodynamic equations and the equa- 
tions of radiation transport. The system is sep- 
arated applying operator splitting, and the par- 



tial equations are then discretized explicitly or im- 
plicitly according to stability considerations. The 
dust size distribution is binned into N mass in- 
tervals and the hydrodynamic equations for mass 
and momentum conservation are computed for the 
gas component and for each binned dust compo- 
nent (grain size) simultaneously. The equations of 
mass conservation for gas {g) and dust particles 
(Qk) can be written: 

|+V.(,v)=0 



dok 
dt 



N 



■ V • (gfcv) = -Qk ^(ttifc + Pik) 



rrii 



^ N N 

+l:^^i(x^jd^jk + liijgijk) — — ■ rrik 

Here, dust coagulation and shattering has been 
included in the equation of continuity for the dust 

particles as a source/sink term. 

The equations for momentum conservation in 
cylindrical coordinates for gas and dust grains are: 

d{QV^) dp i9$ 

^ + V-(,..v) = ---,- 

N 

+ e^Ik{vk,z - Vz) - (V • Qvisc)2 

fe=l 

d{QVr) , ^ , ^ dp 9$ 
^+V-(,..v) = ---,- 



N 



X ^ ^rh 

+ Q2_^Ik{Vk,r - Vr) + Q— - (V • Qvisc) 



k=l 



9{gv4,r) 



dt 



N 



Q'^h{vk,4> -v^) r-{V ■ Qvisc)<;i r 



k=l 

d{QkVk,z) , ^ ( . 'i^k^z 9$ 

— — \- V • {gkVk,z^k) = Qk— Qk^ 

at c oz 

-Qlk{vk,z - Vz) 
-KT^ + V • {QkVk,r^k) = Qk— ^fe-^ 

at c or 



^glk{vk,r - Vr) + Qk-^ 



d{QkVk,4>'r) 



dt 



+ V • {QkVk,4>r^k) = 



-Qh{vk,<j, -V4,) r 
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with the interaction term: 



4 = t^fc — V'c2 + (v-v,)2 (35) 



The advcction part of these equations is solved 
by an exphcit second order scheme. The gas-dust 
interaction terms need an impHcit treatment be- 
cause of the stiffness of the problem. 

The tensor Qvisc which appears in several of the 
above equations denotes the viscous stress tensor 
of the a- viscosity (Shakura & Sunyaev 1973): 

Qvisc = Qve , (36) 
where e is the shear tensor, u is calculated from 



aCsH kQ.T acl/^ , 



(37) 



where the density scale height H has been replaced 
by an expression valid for equilibrium "thin" disks. 
For the calculation described in this study both 
sound speed Cs and angular velocity are ap- 
proximately constant along the surfaces of cylin- 
ders within the equilibrium disk (c.f the theoretical 
discussion in Tassoul 1978). Thus, v varies prin- 
cipally as a function of the radial distance within 
the disk. 

The viscosity as described above is applied to 
the entire grid. Within the accretion disk it modi- 
fies the flow by allowing angular momentum to be 
transfcrcd radially outwards within the disk. The 
parameter a is continuously adjusted according to 
the Toomre stability criterion and is allowed to 
vary within the range of 10"'' to 0.1 (c.f. Yorke 
& Bodenheimer 1999). We define the Toomre pa- 
rameter Qt within the accretion disk for each time 
step by the minimum value of: 



Qt = min 

r < -Rdiak 



=0 ' 



(38) 



where E(r) = / gdz is the disk's surface density. 
If Qt drops below 1.3 (i.e. nonradial instabili- 
ties can be expected to occur), we increase a by 
a small factor (typically 1.002), if necessary to its 
maximum value 0.1. If Qt increases above 1.5 
(the disk becomes stable), a is reduced by a small 
factor (typically 0.999). Generally, Qt 'hovers' at 
either 1.3 or 1.5 and a levels off at values some- 
where between a « 0.01 (for M = 1 Mq) and 
a « 0.08 (for M = 10 Mq). 



Radiation transport is calculated within the 
framework of the grey flux limited diffusion 
approximation (Levermore & Pomraning 1981; 
Yorke & Bodenheimer 1999): 



dt 



(39) 

with flux limiter A^, Rosseland mean opacity 
KkiTd), the radiation constant a, and luminosity 
of the central source L* (treated as an additional 
source term within the volume Vi of the innermost 
grid cell): 

§(yi) — I ■'■/^i innermost cell ^^q^ 
1 otherwise 



coth(^) 



(41) 
(42) 



Here, B^ = Ba(T) is the Planck function. Because 
we are considering grey radiation transfer only, the 
equilibrium temperature of each dust grain is iden- 
tical to the radiation temperature T^. To solve 
equation 39 for Td an implicit ADI procedure is 
used (Douglas & Rachford 1956). 

By contrast, the gas is poorly coupled to the 
radiation field due to its low opacity. The dust 
grains interact with the gas by inelastic collisions 
which contribute a cooling term A to the equation 
for the internal energy density of the gas. The 
dissipation of viscous energy leads to an additional 
gas heating term Qvisc : Vv 

de 

+ V • (ev) = -pV • V - A(^., sk, T, Td) 

+ Qvisc : Vv (44) 



dt 



Assuming an energy transfer of k(Trf — T) per colli- 
sion of a gas molecule with a dust grain the cooling 
function becomes {fio is mean molecular weight of 
the gas): 



He, 8k,T,Td) = 

N 

Q s:^ Qk 



l^orriH ^ ruk 



-crk\ 



I 8kT 
IT Horn H 



k(Td-T) (45) 



The equation of state p{g,T) and the inter- 
nal energy e{g,T) for the gas component assume 
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molecular gas and includes dissociation of the H2 
molecules above « 2000 K (Black & Bodenheimer 
1975). 

Finally, the gravitation potential of the molecu- 
lar cloud is calculated by a solution of the Poisson 
equation, again using ADI: 

N 

A^ = 4ttG {g + J2Qk) . (46) 
fe=i 

4. Initial and boundary conditions 

Table 1: INITIAL CONDITIONS 



Model 


Dust 


[Me] 


[10-12 s-l] 


[yr] 


IMS 


comp. 


1 


1 


8635 


3MS 


comp. 


3 


3 


4985 


SMS 


comp. 


5 


4 


3870 


lOMS 


comp. 


10 


5 


2730 


1MS.PCA 


BPCA 


1 


1 


8635 


3MS.PCA 


BPCA 


3 


3 


4985 


5MS.PCA 


BPCA 


5 


4 


3860 


10MS.PCA 


BPCA 


10 


5 


2730 


IMS.CCA 


BCCA 


1 


1 


8635 


3MS_CCA 


BCCA 


3 


3 


4985 


5MS_CCA 


BCCA 


5 


4 


3860 


lOMS.CCA 


BCCA 


10 


5 


2730 



Note. — Mc: Mass of cloud core; Q: Angulcir velocity; tg: 
Free-fall timescale. 

Wc start the numerical simulations with an 
isothermal, uniformly rotating molecular cloud 
core with a total mass of 1 Mq to 10 Mq, a radius 
of r = 2 X IQi^ cm and a temperature of T = 20 K. 
This gives an initial free-fall timescale tg between 
« 8600 yr and « 2700 yr for these configurations. 
The angular velocities fl range between 10-^2 s-^ 
and 5 x 10-^2 s-^ (see Table 1). We consider cen- 
trally peaked mass distributions £> oc l/(r2 -|- 2;2). 
The total mass contribution of the dust grains is 
set to a fraction of 0.25 x IO-2 of the gas mass (cor- 
responding to the mass contribution of silicates). 

At the outer boundary of the space integration 
domain (\/r2 -|- z"^ = 2 x 10^^ cm) the hydrody- 
namic variables are held constant (no inflow or 
outflow). This corresponds to an assumption that 



no material from the parent cloud will enter into 
the portion undergoing collapse. This does not 
mean that the mass influx rate of dusty material 
onto the new formed disk is suddenly cut off after 
one free-fall time. Instead, it steadily decreases 
as the material in the outer zones, initially slowed 
by pressure gradients in the density-peaked dis- 
tribution, is depleted (see Yorke & Bodenheimer 
1999). Even after three free-fall times, there is an 
appreciable mass influx onto the disk. From other 
studies (e.g. Mizuno et al. 1988) we know that if 
the enshrouding molecular cloud can continuously 
supply small dust grains, the grain size distribu- 
tion will be affected. This effect will be present to 
some extent in the studies discussed here. The is- 
sue of the time dependency of the mass influx is a 
non-trivial one, however, and is beyond the scope 
of the present investigation. 

The particle sizes are distributed according to a 
MRN power law with an exponent of —3.5 which 
can be transformed to a bin mass distribution 
mbm(mfe): 

n{a) oc a-^-^ (47) 
nibin(?n)dlnm oc m ^/^dlnm (48) 

with grain mass m. Figure 4 (upper left panel) 
shows the initial MRN dust mass distribution. On 
a logarithmic scale the dust mass increases with in- 
creasing bin mass. The grain sizes range from 5 nm 
to 5 /im at the beginning. Above 5 /im the particle 
density falls off proportional to mbin~i (this cor- 
responds to n(a) oc a-^). The dynamics of these 
large grains at the upper end of the computed 
size distribution can be followed throughout the 
simulation without being relevant for coagulation 
(provided that the upper end is chosen far enough 
away of the largest particles produced by coagula- 
tion). The integration domain in dust size space 
ranges from 5 nm up to 0.2 mm which spans about 
14 orders of magnitude in grain mass. Whenever 
the dust temperature is low enough to allow an ice 
coating on the grains' surfaces, the sticking propa- 
bilities are modified accordingly. 

The innermost cell of the finest grid contains 
the central protostellar source and requires special 
treatment. Its total luminosity L* can be approx- 
imated by the sum of the core's intrinsic luminos- 
ity and the luminosity due to accretion of material 



10 



onto the core: 

The radius of the central object is held con- 
stant at lOi?0. M*(i) and M*(t) = / M^dt result 
from the calculations. 

Figure 7 (left panel, solid line) displays the net 
specific extinction coefficient for a gas diist mix- 
ture with compact spherical grains. The opaci- 
ties are calculated for the MRN mass distribution 
shown in Figure 4 (upper left panel). Because the 
dust opacities arc calculated using the actual grain 
size distribution they vary as a function of time 
and location during the simulation. 

5. Numerical Simulations 

The following numerical calculations were con- 
ducted with three nested grids of increasing res- 
olution of factor two each. The individual grids 
span 60 x 60 zones. The dust size distribution is 
sampled by = 30 discrete dust species. Table 2 
summarizes selected results of the simulations. 



Table 2: SELECTED RESULTS 













Model 


[Mo] 


[L©] 


Md/M, 


ts/ts 


IMS 


0.77 


3.7 


0.30 


2.5 


3MS 


2.2 


87 


0.36 


2.3 


5MS 


3.7 


194 


0.35 


2.1 


lOMS 


8.2 


1839 


0.22 


1.9 


1MS_PCA 


0.77 


4.0 


0.30 


2.5 


3MS_PCA 


2.3 


59 


0.30 


2.3 


5MS_PCA 


3.7 


270 


0.35 


2.2 


10MS.PCA 


7.9 


1639 


0.27 


2.2 


1MS_CCA 


0.77 


2.9 


0.30 


3.1 


3MS_CCA 


2.2 


80 


0.36 


2.5 


5MS_CCA 


3.8 


188 


0.32 


2.3 


lOMS.CCA 


8.0 


1614 


0.25 


2.2 



Note. — M* and L*: Mass and luminosity of the central 
star; M(j: Mass of disk; tg'- Duration of simulation. 



5.1. Compact spherical dust grains 

The first simulation applies the simple "com- 
pact spherical grain" dust model (section 2.1). 



The collapse of the rotating molecular cloud core 
is followed for about 10** yr (about two initial free- 
fall times). Figure 3 displays an evolutionary se- 
quence of the gas density and the gas velocity and 
Figure 4 the corresponding total grain mass spec- 
trum. 

As evident in the lower right panel of Figure 3 
two accretion shock fronts have developed around 
the protostellar disk. The central mass and core 
luminosity have attained values 2.2 Mq and 87 Lq, 
respectively. As shown in Figure 4, coagulation 
removes the small dust particles effectively from 
the size distribution. At the high mass end dust 
grains of ~ 50 /xm are grown by coagulation. How- 
ever, larger particles do not form during the simu- 
lation. The integral size distribution between 5 nm 
and 5 /xm varies as n{a) oc a^'^'^. 

Figure 5 (left panel) shows the dust mass spec- 
trum at selected heights above the protostellar ac- 
cretion disk for a cylindrical distance of 30 AU. At 
the disk midplane (right panel) coagulation is very 
strong. Large grains are produced at the cost of 
the small size end of the particle spectrum. At disk 
radii larger than about 200 AU the effect of coag- 
ulation becomes less important. Only the small 
grains are removed from the size spectrum by co- 
agulation. However, in the accretion shock just 
above the disk (Fig. 5, left panel) the large grains 
are depleted relative to the low mass dust parti- 
cles. 

The velocities of selected dust grains through 
the accretion shock at r = 30 AU (Fig. 6, left 
panel) show that dust grains of radii ~ 1 /zm 
and above are coupled only loosely to the gas so 
that significant relative velocities of several kms^^ 
between these grains and the smaller ones are 
achieved. This implies that coagulation is inhib- 
ited and shattering might occur (the threshold ve- 
locity is 2.7kms~^). The densities are so low, 
however, that the shattering time scale is 10* yr. 
Thus, the depletion of the large grains must be at- 
tributed to the fact that they pass quickly through 
this zone, whereas the smaller grains slow down 
significantly. This size dependent gas-grain drift 
lowers the dust to gas mass ratio by more than 
a factor of 2 between the two accretion shocks 
(shown for the BPCA fractals, see Fig. 9, right 
panel). 

Note the differing velocities above the outer ac- 
cretion shock {z « 170 AU), caused by the size de- 
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Fig. 3. — Evolution of density and velocity using the compact sphere dust model 3MS. Density contour lines 
are separated by Alog q = 0.5. From left to right and top to bottom the following times are shown: yr, 
5100 yr, 10300 yr and 11400 yr. In the lower right frame the locations of the accretion shocks are given 
by thick dark lines on the left half, where the density contour lines (shown on the right half only) begin 
to tightly bunch. The inner accretion shock encompasses the equilibrium accretion disk, where the {vriVz) 
components of velocity are negligible. The outer accretion shock is less apparant, because the bunching of a 
mere two contour lines — corresponding to a jump in density by a factor of -^10 — is less conspicuous. 
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Fig. 5. — Model SMS, 11400 yr. Left: Dust mass distribution at selected heights above the disk's midplane 
for r = 30 AU. Right: Dust mass distribution at selected radial distances for z = AU. 
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Fig. 6. — Model SMS, 11400 yr. Left: Vertical components of velocity of selected dust components at 
r = 30 AU. Right: Radial velocities of selected dust components in the equatorial plane {z = AU). 




Fig. 7. — Model SMS, 11400 yr. Left: The specific extinction coefficient of the dust in the accretion shock 
{dotted line) and in the equatorial plane of the accretion disk {dashed line) is compared with an MRN mass 
distribution {solid line). Right: Overall mass spectrum at the end of a comparison simulation whereby the 
contribution of systematic velocities to the relative motion of the dust grains is neglected. 
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pendence of the grain opacity and thus by differen- 
tial radiative acceleration. Without radiative ac- 
celeration the grains should fall towards the equa- 
torial plane faster than the gas, because they are 
not pressure supported. Outside the outer accre- 
tion shock the grains' mass distribution is still well 
approximated by the initial MRN mass distribu- 
tion (Fig. 5, left panel). In outer regions of the 
cloud at the equatorial plane the infalling mate- 
rial is shielded from the central star by the disk. 
Hence, the grains do indeed accrete with higher 
velocities than the gas component (Fig. 6, right 
panel). 

To quantify the effect of a modified dust size 
spectrum on the optical properties of the proto- 
stellar matter in the accretion disk. Figure 7 (left 
panel) compares the net specific extinction coef- 
ficient for several locations in the accretion disk. 
Whereas the depletion of the large grains behind 
the accretion shock does introduce minor modi- 
fications to the extinction coefficient, coagulation 
at the disk's midplane causes an opacity reduction 
of more than an order of magnitude for the near 
infrared to UV extinction. From A = 0.1cm to 
0.1mm the extinction coefficient is increased by 
approximately the same amount. This behaviour 
is indicative of the migration of the peak of the 
grain mass distribution to higher masses due to 
coagulation. 

To drive the coagulation beyond the initial up- 
per grain size limit of 5 /im, systematic relative 
velocities are needed, such as those which result 
from differential radiative acceleration, relaxation 
behind the accretion shock and gravitative sed- 
imentation. Figure 7 (right panel) displays the 
resulting dust mass distribution, when these con- 
tributions to the relative motions are neglected. 
Obviously, turbulence and Brownian motion are 
sufficient to remove the small grains, but they are 
not able to build up /xm-sized particles quickly. 
Whereas for cloud clump masses of 1 M0 the dif- 
ferential radiative acceleration of grains is negli- 
gible, it eventually becomes the most important 
mechanism for creating velocity differences be- 
tween grains outside the accretion shock fronts 
when larger clump masses are considered. 



5.2. Influence of the sticking properties of 

the grains 

How is coagulation affected by the sticking 
strength of the dust particles? As pointed out in 
the introduction, this material property is not yet 
well understood. Reports on experimental studies 
indicate larger values for the critical sticking ve- 
locity than theoretical models predict (Poppe & 
Blum 1997). Ices on the grain surfaces play an 
important role (Supulver et al. 1997). To investi- 
gate this effect we conducted comparison simula- 
tions with different sticking strengths. First, the 
critical sticking velocities were reduced to the con- 
servative theoretical values (without the factor 10 
correction for the experimental results, see section 
2.1.2). 

Figure 8 (left panel) displays the total mass 
distribution. Obviously, grain growth is reduced; 
the maximum grain mass is smaller by a factor of 
about 10. When the critical sticking velocity is set 
to infinity, i.e. the grains stick at every encounter, 
grains up to the actual bin limit with radii of about 
0.2 mm are grown. This demonstrates that the 
material parameters play an extremely important 
role in defining the upper mass limit up to which 
the dust grains are able to grow during the forma- 
tion of an accretion disk. 

5.3. Fractal BPCA particles 

Our next approximation treats the dust grains 
as fractal coagulated particles (BPCA particles). 
The same initial and boundary conditions as in 
the previous calculation are used. At the end of 
the simulation the mass of the c;entral object is 
2.3 M0 with a luminosity of 59 Lj . The overall 
evolution is qualitatively similar to the calculation 
with compact dust grains (Fig. 3 and Fig. 4). 
Figure 9 shows the mass density and velocity of 
the gas component (left panel) and the dust to 
gas mass ratio (right panel). Again, an accretion 
shock has developed. In this accretion shock the 
dust to gas mass ratio is reduced by a factor of 
about 3 compared to the initial value due to size 
dependent advection. 

The specific cross section of large BPCA coagu- 
lates is about a factor of 5 larger than for compact 
spheres of the same mass (section 2.2). Only parti- 
cles with radu of about lOfim and larger decouple 
from the motion of the gas flow in the accretion 
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Fig. 8. — Overall mass spectrum of the grains {solid lines) compared to the initial MRN mass distribution 
{dotted lines). Left: The critical sticking velocity is reduced by a factor of 10 compared to model SMS. Right: 
For comparison, the critical sticking velocity is set to infinity. 
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Fig. 9. — Model 3MS_PCA, 11400 yr. Left: Mass density and velocity of the gas component. Symbols and 
contour lines are as in Fig. 3. Right: Variation of the dust to gas mass ratio in the accretion disk. Qi/ q)o 
is the initial value. The outer accretion shock (bunching of ^2 contour lines in the left frame; surface of low 
dust to gas mass ratio in the right frame) &i z k, ±(125 — 150) AU is readily discernible. 
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Fig. 10. Model 3MS_PCA, 11400 yr. Left: Vertical components of velocity of selected dust components at 
r = 30 AU. Right: Overall mass spectrum of the coagulated dust grains {solid line) compared to the initial 
MRN mass distribution {dotted line). 
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Fig. 11. — Model 3MS_CCA, 12600 yr. Left: Gas density and velocity structure. Right: Overall mass 
spectrum of the coagulated dust grains {solid line) compared to the initial MRN mass distribution {dotted 
line). 
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Fig. 12. — Model 3MS_CCA, 12600 yr. Left: Vertical components of velocity of selected dust particles at 
r = 30 AU. Right: The net specific extinction coefficient of unmodified dust {solid line) compared to the 
numerical results for the equatorial midplane at r = 300 AU {dotted line) and r = 30 AU {dashed line). 




Fig. 13.— Model 3MS_CCA_S10, 12700 yr. Left: Mass density and velocity of dust particles with rex = 5 /um. 
Right: Overall mass spectrum of the coagulated dust grains {solid line) compared to the initial MRN mass 
distribution {dotted line). 
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shock (Fig. 10, left panel; compare to Fig. 6). 
Coagulation in the equatorial plane proceeds at a 
rate similar to the simple compact spherical dust 
model (Fig. 10, right panel). As in the case of 
spherical grains no dust particles larger than sev- 
eral 10 /zm are grown by coagulation. This can be 
attributed again to a finite critical sticking veloc- 
ity (see section 5.2). 

5.4. Fractal BCCA particles 

Finally, BCCA grains are the most fluffy parti- 
cles dealt with in these simulations. Figure 11 (left 
panel) shows the density and the velocity of the 
gas component in Model 3MS_CCA. The overall 
distribution is similar to the previously discussed 
models (3MS and 3MS_PCA). However, grain co- 
agulation is enormously strong. Grains as large 
as 0.2 mm (at the limiting end of the binning) are 
grown (Fig. 11, right panel). Almost all the grain 
mass resides in the most massive bin. Because of 
the fluffy structure of the BCCA grains the dust 
remains well coupled to the gas, even in the ac- 
cretion shock (Fig. 12, left panel). Thus, relative 
velocities between the dust particles remain very 
low. 

The optical properties of the gas-dust mixture 
at several locations in the accretion disk are plot- 
ted in Figure 12 (right panel). In the equatorial 
plane the net specific extinction coefficient k°^^ is 
lowered by more than an order of magnitude from 
the near infrared to UV wavelengths. Prom 1 mm 
to 100 /xm the extinction is enhanced. In the ac- 
cretion shock only minor modifications can be as- 
certained. 

A theoretical shattering model differing some- 
what from the one discussed above has been de- 
veloped by Dominik & Tielens (1997). In order 
to test its effect on our simulations we used this 
model to compute the evolution of BCCA grains. 
According to Dominik fc Tielens (1997) the shat- 
tering velocity is proportional to the critical stick- 
ing velocity. It also depends on the number of con- 
tact points between the two colliding dust grains. 
In our ignorance of this quantity we assume that 
two dust grains always have 10 contact points. 
Figure 13 (left panel) displays the density and ve- 
locity of dust grains with reduced radii r^x = 5 /xm. 
In the accretion shock near to the axis of rotation 
these dust particles are destroyed by shattering 
encounters. As can be seen in the total dust mass 



distribution (Fig. 13, right panel), the largest dust 
grains with radii of 0.2 mm are not formed. We at- 
tribute this to frequent shattering encounters; the 
sticking properties are identical to those used in 
model 3MS_CCA. 

5.5. Synthetic emission maps and spectra 

In order to compare our numerical results to ob- 
servations of young protostellar accretion disks we 
have produced emission maps and spectra calcu- 
lated with a 'numerical telescope', which includes 
the contribution of scattered light (Yorke & Bo- 
denheimer 1999). Figure 14 (left panel) shows 
the dust continuum emission at 3.6 /im at the final 
stage of the collapse of the "compact sphere" dust 
model (3MS). A dark bar across the midplane of 
the accretion disk marks a region of high extinc- 
tion. Above and below the disk scattered light 
reveals the presence of the central protostellar ra- 
diation source. The spectral energy distribution 
(SED) shown in Figure 15 (left panel) displays well 
known features of young dusty protostellar cores 
(e.g. Sonnhaltcr, Preibisch & Yorke 1995): No di- 
rect radiation from the protostar (at an angle of 
85°), a silicate absorption feature at A ~ 10 /xm, 
and a dust emission temperature of about 100 K. 
For comparison, the corresponding SED was recal- 
culated for dust with an MRN mass distribution 
using the same density distribution. No obvious 
differences to the results using the more detailed 
coagulated dust model are detectable. This can 
be attributed to the fact that the strongly coag- 
ulated dust particles are embedded deeply within 
the accretion disk, whereas the outer layers con- 
tain only slightly modified grains. However, when 
all dust grains are assumed to have a radius of 
a = 22 /xm (which corresponds to the maximum 
grain size formed in model 3MS), the calculated 
SED is markedly different from that resulting from 
dust with an MRN mass distribution. The silicate 
absorption feature is totally absent and more near 
infrared emission reaches the observer. 

For comparison, an emission map using the un- 
modified MRN dust distribution has been calcu- 
lated and displayed in Figure 14 (right panel). The 
differences are not overwhelming, but some gen- 
eral tendencies can be ascertained: For unmodi- 
fied dust the dark absorption bar across the disk is 
somewhat more prominant (especially towards the 
edges of the disk) and more scattered light above 
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Fig. 14. — Dust continuum emission of model SMS, 11400 yr at an angle of inclination of ^ = 85°. Left: 
Continuum emission map at A = 3.6 /xm using coagulated dust as calculated. Right: Comparison emission 
map using unmodified dust. 




Fig. 15. — Emission spectra resulting from simulations leading to coagulated dust {solid line), with un- 
modified dust (dotted line), or assuming monodisperse dust {dashed line) with a = 22 jiva (model SMS) and 
Tex = 0.2 mm (3MS_CCA). The rotation axis of the disk is inclined at an angle 6 = 85° with respect to the 
line of sight. Left: Model SMS, 11400 yr. Right: Model SMS.CCA, 12600 yr. 
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Fig. 16. — Dust emission of model 3MS_CCA, 12600 yr. The rotation axis is inclined at an angle 6 = 85° 
with respect to the hnc of sight. Left: Continuum emission map at A = 3.6 fim using coagulated dust. Right: 
Emission map assuming unmodified dust. 



and below the disk is visible. Because almost no 
coagulation occurred in these outer disk regions 
over the time period investigated, these differences 
result primarily from the differential advection of 

dust grains. 

A similar tendency can be seen for the frac- 
tal grains. In Figure 16 (left panel) the con- 
tinuum emission for simulation 3MS_CCA is dis- 
played. In contrast to model 3MS the dark absorp- 
tion bar across the disk is far more transparent 
at A = 3.6 /xm when dust coagulation is permit- 
ted. When compared to the emission map result- 
ing from uncoagulated dust with an MRN mass 
distribution (Fig. 16, right panel), it becomes ap- 
parant that the coagulation process has enabled 
the disk to become rather transparent. The over- 
all disk features are in general similar to those re- 
sulting from the simple spherical grain model. 

The SED shown in Figure 15 (right panel) is 
similar to the SED of model 3MS (left panel). The 
SED generated using the coagulated dust from the 
3MS_CCA simulation {solid line, right panel) dis- 
plays some differences with respect to the SED us- 
ing non-coagulated dust: There is a slight shift in 
the emission peak to shorter wavelengths, lower far 
infrared fluxes, enhanced mid-infrared emission, 
and a less prominent silicate absorption feature. 



6. Discussion and Conclusions 

We have shown that dust coagulation pro- 
ceeds at an early phase during the formation of 
a protostellar accretion disk. Small grains with 
0.1 /Km are removed from the mass spectrum 
quickly and effectively in the midplane of the ac- 
cretion disk within ^10^ yr. Here, large particles 
with sizes of several 10 /lui can be produced by co- 
agulation. The maximum grain size which can be 
quickly produced by coagulation during the col- 
lapse and initial accretion of material onto the disk 
depends crucially on the assumed sticking strength 
of the dust particles. In this respect the process 
of ice sublimation plays an important role: When 
the grain surface is coated with material which 
enhances the grain to grain adhesion, the degree 
of coagulation can be significantly increased. 

In the accretion shock front relative velocities 
of several kms"^ arc achieved due to the size de- 
pendent coupling to the gas. Compact spherical 
grains decouple at higher gas densities (and thus 
earlier during the evolution) than fractal dust co- 
agulates. Particle shattering of compact spheri- 
cal grains was not critically important during the 
evolution of the intermediate mass protostars con- 
sidered here. We infer, however, by appropriate 
scaling of masses and luminosities (see Suttner et 
al. 1999), that shattering should be important for 
high mass protostars. For the high mass case ra- 
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diativc acceleration will become increasingly more 
effective in causing a size-dependent spread of dust 
drift velocities. Assuming BCCA grains which 
break apart at even relatively small collision ener- 
gies (as in the model of Dominik & Tielens 1997), 
particle shattering gains some importance for the 
lower cloud masses considered here. Within the 
accretion shock grains are shattered and the max- 
imum grain size is limited to several 10 /zm. How- 
ever, the amoimt of very small debris particles thus 
produced is negligible in the total grain mass spec- 
trum. 

Gas dust drift leads to depletion of dust in the 
immediate vicinity of the accretion disk every- 
where except in the equatorial regions. In par- 
ticular, the gas to dust mass ratio can be low- 
ered by a factor of 2 to 4 within the accretion 
shock. Whereas for a cloud clump mass of 1 M0 
radiative acceleration of dust grains is negligible, 
for clump masses ^ 3 Mq radiative acceleration of 
dust grains becomes increasingly important. De- 
pending on how well the radiation field of the cen- 
tral source is shielded by the disk, the infall of 
dust particles can be hindered in the polar regions, 
whereas in the equatorial regions the dust moves 
radially inwards faster than the gas. 

The optical and physical properties of grains 
are strongly affected by coagulation. The specific 
extinction coefficient in the visual to UV can be 
lowered by more than an order of magnitude in 
the equatorial plane due to coagulation. The grain 
temperature in the midplane and the grains' ca- 
pacity for "freezing out" molecules is correspond- 
ingly affected. Although the local variations of the 
optical coefficients arc large, the only significant 
effect to observational properties is a reduction of 
the near infrared dust opacity in the wavelength 
range between 1 and 100 /im. which is most promi- 
nent for "robust" BCCA particles. Polarization of 
starlight should supply an additional appropriate 
observational tool to determine the degree of co- 
agulation. 

The differences of the global characteristics of 
the simulations using the simple approximation of 
compact spherical grains, BPCA dust, and BCCA 
dust are not as dramatic as may have been naively 
expected. For all three cases the hydrodynamical 
structure (in particular, gas density and velocity) 
is strikingly similar. Thus, we feel justified in us- 
ing our rather crude dust models to perform hy- 



drodynamic simulations of low and intermediate 
mass collapsing clouds and subsequently assume 
more sophisticated detailed dust models to gener- 
ate emission maps, polarization maps, and SEDs. 

Finally, we note that D'Alessio et al. (1999) 
find that synthetic 1 /zm images of accretion disks 
around low mass stars appear to have too large 
geometrical thicknesses to be consistent with ob- 
servation, under the assumption that dust is well 
mixed with the gas. Our study shows that the 
issue might be resolved by taking into proper ac- 
count the differential advection of dust grains. 
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APPENDIX 



Here, we give an analytic expression for the kernel 7(rn, m ,m , Sv{m , m )) in the shattering equation of 
section 2.3. As stated in section 3.1, the integral of the shattering equation is discretized by summing over 
the dust mass space ranging from mi to mjv (or, equivalently, from ai to cat). Thus, 7 transforms to the dis- 
crete function gijk, from which the debris distribution Gk{mi, mj,6vij) = gijkmk/ {mi+mj) can be separated. 



For VcTit < Svij < Vca.t and wrrij e [m^ , m^]: 

Gk {rrii , ruj , Svu ) = — {1 - w) (0-1) 

rrij + m,j 

For Vcrit < Svij < Vcat and Uk < a 

max* 

Gk {mi , rUj , Svij ) = ^ w /mrn (0-2) 

For 6vij > Vce,t and Uk < amax: 



rrii + nij 



Gk {m.i,mj,Svij) = /mrn (0-3) 



Otherwise, Gkimi^nij , Svij) =0. 

We have used the following assumptions and definitions: 



rui > rrij 



m, = {mk+mk-i)/2 = m 



w = 



k - ""k+i 

16/9 

13 



a 



3.64 km s~i 
I'crit = 2.7 km 

Vcat = max^^Wcrit, 1.13 km s~-^ (mj/TOj) 

0.28 (w mj/pbuik)^^^ if Svij < Vcat 

0.20 Ui Vcat /SVij if SVij > Vcat 



n9/16 



The formulae were adapted from the work of Jones et al. (1996). Here, w denotes the ejected crater mass 
in units of the projectile mass, Wcat the critical velocity for the onset of total disruption of the target and 
Omax the radius of the largest debris particle. The debris particles are redistributed according to a MRN size 
distribution /mrn between amin = ai and Omax (i-e. between mmin = ttii and Wmax): 

_ mf/^Amk , . 

■/MRN - ax -5/6. — (0-4) 
E*=i Arm 
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